#######################################################
#######################################################
#######################################################
### CREATED BY JONATHAN KING AND JESSICA SCHOENHERR ###
##### REPLICATION DATA FOR "A MATTER OF OPINION?" #####
#######################################################
#######################################################
#######################################################

library(readxl)
library(tidyverse)
library(ggeffects)
library(stargazer)
library(poliscidata)
library(ggplot2)
library(psych)
library(broom)
library(emmeans)

dataAbortion <- read.csv("AbortionData20210806.csv")

####################
### DATA SORTING ###
####################

# 1 - Abortion, liberal, male
# 2 - Abortion, conservative, male
# 3 - Abortion, liberal, female
# 4 - Abortion, conservative, female
# 5 - Abortion, control

dataAbortion$treatmentGroup <- as.factor(dataAbortion$treatmentGroup)
dataAbortion$decisionTerm <- as.numeric(dataAbortion$decisionTherm)
dataAbortion$partisanship <- as.factor(dataAbortion$partisanship)
dataAbortion$female <- as.factor(dataAbortion$female)

dataAbortion <- within(dataAbortion, treatmentGroup <- relevel(treatmentGroup, ref = 5))

dataAbortion$partisanship <- sjmisc::rec(dataAbortion$partisanship, rec = "1=Democrat;2=Independent;3=Republican", as.num = FALSE)

dataAbortion$scTherm <- as.numeric(dataAbortion$scTherm)
dataAbortion$decisionTherm <- as.numeric(dataAbortion$decisionTherm)

####################
### DESCRIPTIVES ###
####################

dataAbortion %>% group_by(treatmentGroup, partisanship, female) %>% summarize(mean = mean(decisionTherm)) %>% print(n = 30)

dataAbortion$count <- 1

dataAbortion %>% group_by(treatmentGroup, partisanship, female) %>% count(count) %>% print(n = 30)

dataAbortion %>% group_by(partisanship) %>% summarize(mean = mean(decisionTherm))
summary(aov(decisionTherm ~ partisanship, data = dataAbortion))
TukeyHSD(aov(decisionTherm ~ partisanship, data = dataAbortion))
abortPart <- lm(decisionTherm ~ partisanship, data = dataAbortion)
car::Anova(abortPart)

dataAbortion %>% group_by(female) %>% summarize(mean = mean(decisionTherm))
summary(aov(decisionTherm~female, data = dataAbortion))
TukeyHSD(aov(decisionTherm~female, data = dataAbortion))
abortFem <- lm(decisionTherm ~ female, data = dataAbortion)
car::Anova(abortFem)

dataAbortion %>% group_by(treatmentGroup) %>% summarize(mean = mean(decisionTherm))

################
### BASELINE ###
################

abortionTreatment <- lm(decisionTherm ~ treatmentGroup, 
						data = dataAbortion)
summary(abortionTreatment)
nobs(abortionTreatment)

stargazer(abortionTreatment, type = "text")

car::Anova(abortionTreatment)

#abortionCoef <- tidy(abortionTreatment)
#abortionCoef$upperCI <- (abortionCoef$std.error * 1.96) + abortionCoef$estimate
#abortionCoef$lowerCI <- (abortionCoef$std.error * -1.96) + abortionCoef$estimate

abortionTreatmentPreds <- emmeans::emmeans(abortionTreatment, specs="treatmentGroup")
pairs(abortionTreatmentPreds)
abortionTreatmentPreds <- as.data.frame(abortionTreatmentPreds)

abortionTreatmentPlot <- ggplot(abortionTreatmentPreds) +
	geom_bar( aes(x = treatmentGroup, y = emmean), stat="identity", fill = "grey75") +
	geom_errorbar( aes(x = treatmentGroup, ymin = lower.CL, ymax = upper.CL), width=0.1) +
	theme_bw() +
	xlab("\nOpinion Writer Treatment Group")  +
	ylab("Paricipants' Warmth Toward Decision\n") +
	ggtitle("Abortion Vignette - Direct Treatment Effects") +
	theme(plot.title = element_text(hjust = 0.5)) +
	scale_x_discrete(labels = c("Control", "Liberal\nMale\nJustice", "Conservative\nMale\nJustice", "Liberal\nFemale\nJustice", "Conservative\nFemale\nJustice")) +
	scale_y_continuous(limits = c(0, 80), breaks = seq(0, 80, 10))
abortionTreatmentPlot

#abortionTreatmentCoef <- ggplot(abortionCoef[2:5,], aes(y = estimate, x = term, ymin = lowerCI, ymax = upperCI)) + 
#	geom_pointrange(fatten = 2, size = 1) + 
#	theme_bw() +
#	xlab("\nOpinion Writer Treatment Group") +
#	ylab("Difference in Support Compared to Control\n") +
#	ggtitle("Abortion Vignette Coefficient Plot - Baseline Model") +
#	theme(plot.title = element_text(hjust = 0.5)) +
#	scale_x_discrete(labels = c("Liberal\nMale\nJustice", "Conservative\nMale\nJustice", "Liberal\nFemale\nJustice", "Conservative\nFemale\nJustice")) + 
#	scale_y_continuous(limits = c(-20, 10), breaks = seq(-20, 10, 5)) +
#	theme(legend.position = "none") +
#	geom_hline(yintercept = 0, size = 1, linetype = "dashed")
#abortionTreatmentCoef

################################
### CONTROL FOR PARTISANSHIP ###
################################

abortionPartisanship <- lm(decisionTherm ~ treatmentGroup * partisanship,
						data = dataAbortion)
summary(abortionPartisanship)
nobs(abortionPartisanship)

car::Anova(abortionPartisanship)

abortionPartisanshipTestDem <- data.frame(treatmentGroup = c(factor(1), factor(2), factor(3), factor(4), factor(5)), partisanship = "Democrat")
predict(abortionPartisanship, newdata = abortionPartisanshipTestDem)

abortionPartisanshipTestInd <- data.frame(treatmentGroup = c(factor(1), factor(2), factor(3), factor(4), factor(5)), partisanship = "Independent")
predict(abortionPartisanship, newdata = abortionPartisanshipTestInd)

abortionPartisanshipTestRep <- data.frame(treatmentGroup = c(factor(1), factor(2), factor(3), factor(4), factor(5)), partisanship = "Republican")
predict(abortionPartisanship, newdata = abortionPartisanshipTestRep)

abortionPartisanshipPreds <- emmeans::emmeans(abortionPartisanship, specs = "partisanship", by = "treatmentGroup")
pairs(abortionPartisanshipPreds)

abortionPartisanshipPredsDF <- as.data.frame(abortionPartisanshipPreds)

aPartisanPlot <- ggplot(abortionPartisanshipPredsDF, aes(y = emmean, x = treatmentGroup, color = partisanship, shape = partisanship)) +
	geom_point(position = position_dodge(0.3)) + 
	geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0, position = position_dodge(0.3)) +
	theme_bw() +
	xlab("\nOpinion Writer Treatment Group") +
	ylab("Predicted Value - Decision Thermometer\n") +
	ggtitle("Abortion Vignettes by Participant Partisanship") +
	theme(plot.title = element_text(hjust = 0.5)) +
	scale_x_discrete(labels = c("Control", "Liberal\nMale\nJustice", "Conservative\nMale\nJustice", "Liberal\nFemale\nJustice", "Conservative\nFemale\nJustice")) + 
	scale_y_continuous(limits = c(20, 100), breaks = seq(20, 100, 10)) +
	theme(legend.position = "bottom", legend.title = element_blank()) +
	scale_color_manual(values = c("black", "grey65", "grey50"), labels = c("Democrat", "Independent", "Republican")) + 
	scale_shape_manual(values = c(15, 17, 19), labels = c("Democrat", "Independent", "Republican"))
aPartisanPlot

##########################
### CONTROL FOR GENDER ###
##########################

abortionGender <- lm(decisionTherm ~ treatmentGroup * female, 
				data = dataAbortion)
summary(abortionGender)
nobs(abortionGender)

car::Anova(abortionGender)

abortionGenderTestMale <- data.frame(treatmentGroup = c(factor(1), factor(2), factor(3), factor(4), factor(5)), female = factor(0))
predict(abortionGender, newdata = abortionGenderTestMale)

abortionGenderTestFemale <- data.frame(treatmentGroup = c(factor(1), factor(2), factor(3), factor(4), factor(5)), female = factor(1))
predict(abortionGender, newdata = abortionGenderTestFemale)

abortionGenderPreds <- emmeans::emmeans(abortionGender, specs = "female", by = "treatmentGroup")
pairs(abortionGenderPreds)

abortionGenderPredsDF <- as.data.frame(abortionGenderPreds)

abortionGenderPreds2 <- emmeans::emmeans(abortionGender, specs = "treatmentGroup", by = "female")
pairs(abortionGenderPreds2)

aGenderPlot <- ggplot(abortionGenderPredsDF, aes(y = emmean, x = treatmentGroup, color = female, shape = female)) +
	geom_point(position = position_dodge(0.3)) + 
	geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0, position = position_dodge(0.3)) +
	theme_bw() +
	xlab("\nOpinion Writer Treatment Group") +
	ylab("Predicted Value - Decision Thermometer\n") +
	ggtitle("Abortion Vignettes by Participant Gender") +
	theme(plot.title = element_text(hjust = 0.5)) +
	scale_x_discrete(labels = c("Control", "Liberal\nMale\nJustice", "Conservative\nMale\nJustice", "Liberal\nFemale\nJustice", "Conservative\nFemale\nJustice")) + 
	scale_y_continuous(limits = c(20, 100), breaks = seq(20, 100, 10)) +
	theme(legend.position = "bottom", legend.title = element_blank()) +
	scale_color_manual(values = c("black", "grey50"), labels = c("Male", "Female")) +
	scale_shape_manual(values = c(15, 19), labels = c("Male", "Female"))
aGenderPlot

###############################
### GENDER AND PARTISANSHIP ###
###############################

modelAbortion <- lm(decisionTherm ~ treatmentGroup
						+ female
						+ partisanship
						+ female * partisanship
						+ treatmentGroup * female * partisanship
						+ treatmentGroup * female
						+ treatmentGroup * partisanship, 
						data = dataAbortion)
summary(modelAbortion)
nobs(modelAbortion)

stargazer(modelAbortion, type = "text")

car::Anova(modelAbortion)

abortionMaleDem <- data.frame(treatmentGroup = c(factor(1), factor(2), factor(3), factor(4), factor(5)), female = factor(0), partisanship = "Democrat")
predict(modelAbortion, newdata = abortionMaleDem)

abortionMaleInd <- data.frame(treatmentGroup = c(factor(1), factor(2), factor(3), factor(4), factor(5)), female = factor(0), partisanship = "Independent")
predict(modelAbortion, newdata = abortionMaleInd)

abortionMaleRep <- data.frame(treatmentGroup = c(factor(1), factor(2), factor(3), factor(4), factor(5)), female = factor(0), partisanship = "Republican")
predict(modelAbortion, newdata = abortionMaleRep)

abortionFemaleDem <- data.frame(treatmentGroup = c(factor(1), factor(2), factor(3), factor(4), factor(5)), female = factor(1), partisanship = "Democrat")
predict(modelAbortion, newdata = abortionFemaleDem)

abortionFemaleInd <- data.frame(treatmentGroup = c(factor(1), factor(2), factor(3), factor(4), factor(5)), female = factor(1), partisanship = "Republican")
predict(modelAbortion, newdata = abortionFemaleInd)

abortionFemaleRep <- data.frame(treatmentGroup = c(factor(1), factor(2), factor(3), factor(4), factor(5)), female = factor(1), partisanship = "Independent")
predict(modelAbortion, newdata = abortionFemaleRep)

abortionPreds <- emmeans::emmeans(modelAbortion, specs = c("female", "partisanship"), by = "treatmentGroup")
pairs(abortionPreds)

abortionPreds <- emmeans::emmeans(modelAbortion, specs = "treatmentGroup", by = c("female", "partisanship"))
pairs(abortionPreds)
	
abortionPredsDF <- as.data.frame(abortionPreds)

abortionPredsDF$group <- ifelse(abortionPredsDF$partisanship == "Democrat" & abortionPredsDF$female == 0, 1, 
ifelse(abortionPredsDF$partisanship == "Democrat" & abortionPredsDF$female == 1, 2,
ifelse(abortionPredsDF$partisanship == "Independent" & abortionPredsDF$female == 0, 3,
ifelse(abortionPredsDF$partisanship == "Independent" & abortionPredsDF$female == 1, 4,
ifelse(abortionPredsDF$partisanship == "Republican" & abortionPredsDF$female == 0, 5,
ifelse(abortionPredsDF$partisanship == "Republican" & abortionPredsDF$female == 1, 6, 0))))))
abortionPredsDF$group <- as.factor(abortionPredsDF$group)

#library(jtools)
#plot_coefs(modelAbortion)

#abNames <- as_labeller(c(`1`="Male Democrat Participants", `2`="Female Democrat Participants", `3`="Male Independent Participants", `4`="Female Independent Participants", `5`="Male Republican Participants", `6`="Female Republican Participants"))

#abortionPlot1 <- ggplot(abortionPredsDF, aes(y = emmean, x = treatmentGroup, ymin = lower.CL, ymax = upper.CL)) +
#	geom_pointrange(fatten = 2, size = 1) + 
#	geom_point(position = position_dodge(0.3)) + 
#	facet_wrap(~ group, scales = "free", labeller = abNames) + 
#	theme_bw() +
#	xlab("\nOpinion Writer Treatment Group") +
#	ylab("Predicted Value - Decision Thermometer\n") +
#	ggtitle("Abortion Vignette Full Results") +
#	theme(plot.title = element_text(hjust = 0.5)) +
#	scale_x_discrete(labels = c("Control", "Liberal\nMale\nJustice", "Conservative\nMale\nJustice", "Liberal\nFemale\nJustice", "Conservative\nFemale\nJustice")) + 
#	scale_y_continuous(limits = c(10, 90), breaks = seq(10, 90, 20)) +
#	theme(legend.position = "bottom", legend.title = element_blank())
#abortionPlot1

#abortionPreds2 <- emmeans::emmeans(modelAbortion, specs = c("treatmentGroup", "female"), by = "partisanship")
#pairs(abortionPreds2)

#abortionPredsDF$treatmentGroup <- factor(abortionPredsDF$treatmentGroup, levels = c(5, 1, 2, 3, 4))

#abNames2 <- as_labeller(c(`5`="Control", `1`="Liberal Male Justice Treatment", `2`="Cons. Male Justice Treatment", `3`="Liberal Female Justice Treatment", `4`="Cons. Female Justice Treatment"))

#abortionPlot2 <- ggplot(abortionPredsDF, aes(y = emmean, x = partisanship, ymin = lower.CL, ymax = upper.CL, color = female, shape = female)) +
#	geom_pointrange(fatten = 2, size = 1, position = position_dodge(0.4)) + 
#	facet_grid(~treatmentGroup, labeller = abNames2) + 
#	theme_bw() +
#	xlab("\nParticipant Partisanship") +
#	ylab("Predicted Value - Decision Thermometer\n") +
#	ggtitle("Abortion Vignette Expanded Results") +
#	theme(plot.title = element_text(hjust = 0.5)) +
#	scale_x_discrete(labels = c("Dem", "Ind", "Rep")) + 
#	scale_y_continuous(limits = c(10, 90), breaks = seq(10, 90, 10)) +
#	theme(legend.position = "bottom", legend.title = element_blank()) +
#	scale_color_manual(values = c("black", "grey50"), labels = c("Male", "Female")) +
#	scale_shape_manual(values = c(15, 19), labels = c("Male", "Female"))
#abortionPlot2

abNames3 <- as_labeller(c(`Democrat` = "Democrat Participants", `Independent` = "Independent Participants", `Republican` = "Republican Participants"))

#abortionPlot3 <- ggplot(abortionPredsDF, aes(y = emmean, x = treatmentGroup, ymin = lower.CL, ymax = upper.CL, color = female, shape = female)) +
#	geom_pointrange(fatten = 2, size = 1, position = position_dodge(0.4)) + 
#	facet_grid(~partisanship, labeller = abNames3) + 
#	theme_bw() +
#	xlab("\nOpinion Writer Treatment Groups") +
#	ylab("Predicted Value - Decision Thermometer\n") +
#	ggtitle("Abortion Vignette Expanded Results") +
#	theme(plot.title = element_text(hjust = 0.5)) +
#	scale_x_discrete(labels = c("Control", "Liberal\nMale\nJustice", "Conservative\nMale\nJustice", "Liberal\nFemale\nJustice", "Conservative\nFemale\nJustice")) + 
#	scale_y_continuous(limits = c(10, 100), breaks = seq(10, 100, 10)) +
#	theme(legend.position = "bottom", legend.title = element_blank()) +
#	scale_color_manual(values = c("black", "grey50"), labels = c("Male", "Female")) +
#	scale_shape_manual(values = c(15, 19), labels = c("Male", "Female"))
#abortionPlot3

abortionPlot4 <- ggplot(abortionPredsDF, aes(y = emmean, x = treatmentGroup, ymin = lower.CL, ymax = upper.CL, color = female, fill = female)) +
	geom_bar(position = position_dodge(), stat = "identity") +
	facet_wrap(~partisanship) + 
	geom_errorbar(position = position_dodge(0.9), width=0.2, color = "grey25") + 
	theme_bw() +
	xlab("\nOpinion Writer Treatment Groups") +
	ylab("Predicted Value - Decision Thermometer\n") +
	ggtitle("Abortion Vignette Expanded Results") +
	theme(plot.title = element_text(hjust = 0.5)) +
	scale_x_discrete(labels = c("Control", "Liberal\nMale\nJustice", "Conservative\nMale\nJustice", "Liberal\nFemale\nJustice", "Conservative\nFemale\nJustice")) + 
	scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10)) +
	theme(legend.position = "bottom", legend.title = element_blank()) +
	scale_color_manual(values = c("grey75", "grey50"), labels = c("Male", "Female")) +
	scale_fill_manual(values = c("grey75", "grey50"), labels = c("Male", "Female"))
abortionPlot4

#################
### stargazer ###
#################

stargazer(abortionTreatment, abortionGender, abortionPartisanship, modelAbortion, title = "Ordinary Least Squares Regression Results, Decision Thermometer, Abortion Vignettes", align = TRUE)

####################
### SMALL SAMPLE ###
####################

# check to see if there's a difference with the reminder
# results hold if you take the small sample out
largerSampleAbortionData <- dataAbortion %>% filter(smallSample != 1)

modelAbortionTreatmentLarge <- lm(decisionTherm ~ treatmentGroup, 
						data = largerSampleAbortionData)
summary(modelAbortionTreatmentLarge)
nobs(modelAbortionTreatmentLarge)

##############################
### CONTROL FOR LEGITIMACY ###
##############################

modelAbortionLegit <- lm(decisionTherm ~ treatmentGroup
							+ legitScore, 
						data = dataAbortion)
summary(modelAbortionLegit)
nobs(modelAbortionLegit)

modelAbortionLegit2 <- lm(decisionTherm ~ treatmentGroup
							+ legitAverage, 
						data = dataAbortion)
summary(modelAbortionLegit2)
nobs(modelAbortionLegit2)

# legit score is not a significant predictor

##############################
### BIG MODEL FOR APPENDIX ###
##############################

abortionWithControls <- lm(decisionTherm ~ treatmentGroup
						+ female
						+ partisanship
						+ income
						+ education
						+ age
						+ female * partisanship
						+ treatmentGroup * female * partisanship
						+ treatmentGroup * female
						+ treatmentGroup * partisanship, 
						data = dataAbortion)
summary(abortionWithControls)
nobs(abortionWithControls)

abortionBigPreds <- emmeans::emmeans(abortionWithControls, specs = c("female", "partisanship"), by = "treatmentGroup")
pairs(abortionBigPreds)
abortionBigPredsDF <- as.data.frame(abortionBigPreds)

abortionBigPredsDF$treatmentGroup <- factor(abortionBigPredsDF$treatmentGroup, levels = c(5, 1, 2, 3, 4))

names <- as_labeller(c(`5`="Control", `1`="Liberal Male Justice", `2`="Conservative Male Justice", `3`="Liberal Female Justice", `4`="Conservative Female Justice"))

abortionBigPlot <- ggplot(abortionBigPredsDF, aes(y = emmean, x = partisanship, color = female, shape = female)) +
	geom_point(position = position_dodge(0.3)) + 
	geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0, position = position_dodge(0.3)) +
	facet_grid(~treatmentGroup, labeller = names) + 
	theme_bw() +
	xlab("\nParticipant Partisanship") +
	ylab("Predicted Value - Decision Thermometer\n") +
	ggtitle("Abortion Vignette Full Results") +
	theme(plot.title = element_text(hjust = 0.5)) +
	scale_x_discrete(labels = c("Dem", "Ind", "Rep")) + 
	scale_y_continuous(limits = c(00, 100), breaks = seq(00, 100, 10)) +
	theme(legend.position = "bottom", legend.title = element_blank()) +
	scale_color_manual(values = c("black", "grey50"), labels = c("Male", "Female")) +
	scale_shape_manual(values = c(15, 19), labels = c("Male", "Female"))
abortionBigPlot

stargazer(abortionTreatment, abortionGender, abortionPartisanship, modelAbortion, abortionWithControls, title = "Ordinary Least Squares Regression Results, Decision Thermometer, Abortion Vignettes", align = TRUE)

#######################
### ABORTION BACKUP ###
#######################

# try it with the other DV
dataAbortion$agreeDV <- ifelse(dataAbortion$agreeWithDecisionRecode == 1, 1, 0)
dataAbortion <- dataAbortion %>% filter(!is.na(agreeDV))

modelAgreeAbortion <- glm(agreeDV ~ treatmentGroup
						+ female
						+ partisanship
						+ female * partisanship
						+ treatmentGroup * female * partisanship
						+ treatmentGroup * female
						+ treatmentGroup * partisanship, 
						family = binomial(link = logit),
						data = dataAbortion)
summary(modelAgreeAbortion)
nobs(modelAgreeAbortion)

aPred2 <- ggpredict(modelAgreeAbortion, terms = c("treatmentGroup [1, 2, 3, 4, 5]", "female [0,1]", "partisanship [Democrat, Republican]"))

abortionPlot2 <- ggplot(aPred2, aes(y = predicted, x = x, color = group)) +
	geom_point(position = position_dodge(0.3)) + 
	geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, position = position_dodge(0.3)) +
	facet_wrap(~facet) + 
	theme_bw() +
	xlab("\nOpinion Writer Treatment Group") +
	ylab("Predicted Value - Decision Thermometer\n") +
	ggtitle("Abortion Vignette") +
	theme(plot.title = element_text(hjust = 0.5)) +
	scale_x_discrete(labels = c("Liberal\nMale", "Conservative\nMale", "Liberal\nFemale", "Conservative\nFemale", "Control")) + 
	scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, 0.1)) +
	theme(legend.position = "bottom", legend.title = element_blank()) +
	scale_color_manual(values = c("black", "azure4"), labels = c("Male", "Female"))
abortionPlot2

# try it with scTherm as the DV
modelAbortionSC <- lm(scTherm ~ treatmentGroup
						+ female
						+ partisanship
						+ female * partisanship
						+ treatmentGroup * female * partisanship
						+ treatmentGroup * female
						+ treatmentGroup * partisanship, 
						data = dataAbortion)
summary(modelAbortionSC)
nobs(modelAbortionSC)

aPred3 <- ggpredict(modelAbortionSC, terms = c("treatmentGroup [1, 2, 3, 4, 5]", "female [0,1]", "partisanship [Democrat, Republican]", "abortionSupport [1]"))

abortionPlot3 <- ggplot(aPred3, aes(y = predicted, x = x, color = group)) +
	geom_point(position = position_dodge(0.3)) + 
	geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0, position = position_dodge(0.3)) +
	facet_wrap(~facet) + 
	theme_bw() +
	xlab("\nOpinion Writer Treatment Group") +
	ylab("Predicted Value - Decision Thermometer\n") +
	ggtitle("Abortion Vignette") +
	theme(plot.title = element_text(hjust = 0.5)) +
	scale_x_discrete(labels = c("Liberal\nMale", "Conservative\nMale", "Liberal\nFemale", "Conservative\nFemale", "Control")) + 
	scale_y_continuous(limits = c(20, 115), breaks = seq(20, 115, 10)) +
	theme(legend.position = "bottom", legend.title = element_blank()) +
	scale_color_manual(values = c("black", "azure4"), labels = c("Male", "Female"))
abortionPlot3

